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SUMMARY 


A simple clustering transformation is combined with the Thompson, Thames, 
and Mastin (TTM) method of generating computational grids to produce con- 
trolled mesh spacings. For various practical grids, the resulting hybrid 
scheme is easier to apply than the inhomogeneous clustering terms included in 
the TTM method for this purpose. The technique is illustrated in application 
to airfoil problems, and listings of a FORTRAN computer code for this usage 
are included. 


INTRODUCTION 


Thompson, Thames, and Mastin (TTM) have developed a practical way of 
generating two-dimensional computational grids about arbitrary bodies with 
coordinate lines coincident with all boundaries (refs. 1,2). Curvilinear 
nonorthogonal coordinates are generated by conventional numerical solution of 
two nonlinear elliptic partial differential equations. Because a refined grid 
is generally needed in regions where a rapid change in gradient of the flow 
variables occurs, the TTM method includes a provision for concentrating coor- 
dinate lines in any region of the field. However, this clustering feature can 
be either complicated or unsatisfactory, at least for the uninitiated user. 

The purpose of this note is to illustrate a simple alternative procedure for 
clustering about any grid line. 

The authors wish to acknowledge the helpful discussions with L. B. Schiff 
with whom we first investigated the TTM scheme, and the assistance of 
P. Kutler who later programmed an airfoil grid generation scheme (some por- 
tions of which we adapted) to investigate the effectiveness of the P and Q 
terms. The authors also wish to acknowledge the assistance given by 
C. Schulbach in writing computer code for the SC4020 plotter. 


THE THOMPSON-THAMES-MASTIN METHOD 


In its simplest form, the TTM method solves Laplace equations in the 
physical plane 


^ xx + ^ yy ® 

+ n yy 


0 


( 2 ) 



to generate curvilinear coordinates E,(x,y) and n (x,y) , which form a rectangu- 
lar grid in the transformed computational plane. Equations (1) and (2) are 
transformed into the computational plane as 

ax^ - 26^ + Y x nn = 0 (3) 

a % " + Yy nn = 0 (4) 

where a = x n 2 + z/ n 2 , 3 = and Y = x ^ + h ^ • Although equa- 

tions (3) and (4) are mildly nonlinear, it is usually a simple matter to solve 
for x and y on the rectangular £,n computational grid by conventional relaxa- 
tion methods. Then, if the transformation Jacobian at each point is nonzero, 
the locally one-to-one mapping between the computational plane and the physi- 
cal plane yields the solution for £,n in the physical plane as well. 

Simple boundary conditions usually accompany equations (3) and (4) . For 
example, Dirichlet conditions can be picked for the blunt-body problem 
sketched in figure 1, where lines of constant n are chosen to coincide with 
the body and outer shock boundaries and the end points of lines of constant £ 
are arbitrarily spaced along these boundaries. Lines of constant £ terminate 
the blunt body flow region, with end points of lines of constant n distributed 
along these boundaries. This ability to arbitrarily locate boundary points by 
specifying their locations in the x ,y plane is one of the most desirable and 
significant aspects of the TTM method. 

While the Laplace equations automatically generate grids, the nodes are 
often poorly concentrated in regions where a large change of gradient in the 
flow variables occurs. To illustrate this, a grid is constructed about an 
NACA 0012 airfoil using Dirichlet boundary conditions as in figure 2. The 
resulting grid is shown in figure 3(a) , and the lines nearest the airfoil are 
shown in detail in figure 3(b). The computed grid has the deficiency that the 
spacing between the airfoil-wake boundary and the next outward line of con- 
stant n is too large and not uniform, especially near the airfoil trailing 
edge. Clustering of lines of constant q to the airfoil-wake boundary 
( abode in fig. 2) is needed. 

The provision for clustering in the TTM method consists of adding 
inhomogeneous terms to equations (1) and (2), producing Poisson equations 


^ xx + ^ yy (5) 

*xx + n yy = G(£,n) 

which in the computational plane become 

ax^ - 2$x^ + YX nn = -J 2 (Px^ + Qx n ) (7) 

~ 2 + yy nn = -J z (Py^ + QyJ (8) 
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where J = x^y^ - x^y ^ is the transformation Jacobian. Here P and Q are sum- 
mations of exponential terms, as defined in reference 2 (however, the sign of 
these terms in ref. 2 is in error). 

The harmonic functions P and Q can cause clustering about any number of 
nodes or grid lines. The TTM method of clustering is versatile (with respect 
to where clustering takes place in the field) , but exact control of mesh 
spacings can be difficult. Particularly difficult is the imposition of a pre- 
determined and uniform spacing between two adjacent grid lines along their 
entire length. 

To illustrate this last remark, consider the 0012 airfoil application. 
Here terms of Q giving clustering near the airfoil-wake boundary line of con- 
stant q are used. The closeup of the resultant grid, shown in figure 4, is 
improved over the unclustered grid shown in figure 3(b) but is still not suf- 
ficiently uniform around the airfoil. 


COMBINING SIMPLE STRETCHING WITH THE TTM METHOD 


In this section we show that for particular problems in which only grid 
refinement near a boundary is desired, a simple stretching procedure combined 
with the simplest, unclustered form of the TTM method can produce good 
results. This technique is easier to implement than the inhomogeneous form of 
the TTM method. 

In the airfoil problem, the £ spacing can be sufficiently controlled by 
placement in the x,y plane of the nodes along the q min and boundaries, so 

only clustering of lines of constant q is needed. Instead of trying to 
enforce a specific clustering by choice of Q, the lines of constant £ from the 
inhomogeneous form of the TTM method can be accepted and the lines of con- 
stant q can be discarded. Then points on lines of constant n are redistrib- 
uted along each line of constant £ by means of simple monotonic stretching 
functions. This method is detailed below for an airfoil application where an 
exponential stretching function is chosen to give a specified minimum grid 
spacing adjacent to the airfoil-wake boundary. 

Given the grid resulting from the simplest form of the TTM method, the 
arc length 5(n)|^ along each constant £ line is computed from the relations 

S 1 “ 0 (9) 

Sfc = + ~^^ x k ~ X k- 1^ 2 + ^k ~ ^k- 1^ 2 < k < fcjnax (10) 

where k is the index of the ?cth line of constant q. Along each line of con- 
stant £, equations (9) and (10) allow tabulation of xj^ and yfc as functions of 
Sy. The q spacing from the TTM solution can now be discarded. 
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An exponential stretching transformation is introduced such that S(n)|£; 
will be the arc length along the constant £ line with grid-point distribution 
defined by 


s 1 = o 


(ID 

~ S k ■ S k - 1 + 4S «in (1 + 

2 * k S W 

(12) 


where A5’ m £ n is the desired minimum physical spacing along lines of constant E, 

near the body, and e is chosen so that Sb = Sy . Values of xr, and yr. 

max max 

corresponding to the new Sfc are obtained by interpolating the x = x(S0 and 
y = tables generated with equations (9) and (10). These locations 

define the new lines of constant n* A new grid has therefore been defined 
where the lines of constant n are uniformly spaced in the physical plane and 
equal to AS mi - n along the entire airfoil-wake boundary and exponentially 
stretched toward the outer (n = n max ) boundary. 


For each constant E, line, the value of e is determined such that 

Sv = Sb using a Newton-Raphson iterative procedure. With e„ defined as 
''•max K max at- t- n 

the value of e after the wth iteration 


where 


i- v t— , / 

n - 1 

e n = e n-i-F’( Vl ) 

F(e) = S' -~S 

^max ^max 


Combining equations (11), (12), and (14) and reducing yields 


(13) 


(14) 


lAkJ • 7 

F(t) = S ^ [(1 + e ) ""max - 1 _ i] (15) 

K max e 

from which it is seen that 

^min | k - 2 I 

F' (e) = — p- jl + (1 + e)' inax [z(k max - 2) - 1]J . (16) 

A value for Eq on the first constant E, line for which the Newton-Raphson pro- 
cedure has been found to be stable is obtained from 


^/(^max 2 ) 


max 


e 0 = 


.AS . 

1 mm/ 


- 1 


(17) 


For all subsequent constant E, lines, the final solution e for the immediately 
previous line is used for e Q . 
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The above procedure has been applied to the NACA 0012 airfoil with A5 m i n 
as small as 0.00005 chords. An example with AS^-j^ = 0.01 chords is shown in 
figure 5(a) and a closeup of the lines nearest the airfoil is seen in fig- 
ure 5(b). The control of the grid line spacing is evident. In comparing 
figure 5(a) with figure 3(a), one observes that the outermost line of con- 
stant n (in the shape of three sides of a rectangle) has been discarded and 
the next inner line of constant n is taken for the p = n max boundary. This 
was done before application of the simple clustering procedure for the purpose 
of minimizing variations in the lengths of different constant E, lines. The 
appendix contains a listing of the FORTRAN computer code used to generate 
these grids and notes on its use. 


CONCLUDING REMARKS 


A simple stretching transformation has been combined with the simplest 
form of the Thompson- Thames-Mastin grid generation method to produce grids 
with certain types of controlled spacings. The technique has been used to 
generate various grids including those for viscous flow calculations about 
airfoils. 

Ultimately, the clusterings in both E, and q will probably be done by an 
improved version of the Thompson-Thames-Mastin method. Such an algorithm 
should include a procedure for dynamic adjustment of the parameters in the 
P and Q terms during solution of the Poisson equations. For the present, the 
much simpler method is adequate for many cases. 

Ames Research Center 

National Aeronautics and Space Administration 

Moffett Field, Calif. 94035, March 8, 1977 
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APPENDIX 


FORTRAN PROGRAM FOR AIRFOIL GRIDS 


The FORTRAN computer program that combines simple exponential stretching 
with the TTM method to produce two-dimensional grids about airfoils is 
designed to be run on the Control Data Corporation 7600 computer system but 
should prove easily adaptable to other installations. It is approximately 
700 cards in length and consists of a main program and twelve short subrou- 
tines. Of these subroutines, five deal with plotting the resultant grids on 
the Stromberg-Carlson SC4020 plotter and can be easily removed if this feature 
is not desired. Another two of the included decks are proven utility subrou- 
tines for standard mathematical operations. Thus, only the main program and 
five of the subroutines constitute the operative part of the code that creates 
the grids. 


MAIN PROGRAM 


The main program first reads input data cards, which are detailed in a 
later section of this appendix. Subroutine INNER is called to distribute 
points on the airfoil-wake and rear boundaries. Points on the bottom-front- 
top boundary are located by subroutine OUTER. With all boundary points fixed, 
initial conditions for the TTM method are provided by equally spacing grid 
nodes along straight lines of constant £ between corresponding boundary points. 
The TTM method is applied by subroutine RELAX, exponential stretching is 
applied by subroutine CLUSTR, and the finally resultant grid is output using 
a FORTRAN unformatted WRITE statement. 


♦ DECK 

C 

C 

C 

C 

c 

c 

c 

c 

c 


c 

c 

c 

c 

c 

c 


MAIN 

PROGRAM M AIN < INPUT ,TA PE I NPUT, OUTP UT, T APE6 *0UT PUT , TAPE48, TAPES) 
LOGICAL UNIT NUMBERS USED IN THIS PROGRAM* 

5 - READ INPUT DATA CARDS 

6 - WRITE PRINTER OUTPUT 

6 - UNFORMATTED WRITE OF FINAL SOLUTION GRID 

48 - SC4020 PLOTTER DATA 

COMMON/ TCQM/ IT IT LE ( 8 ) 

DIMENSION X8 (140 >,YBU40) ,X( 140, 60) , Y ( 140, 80 ) 

DIMENSION LIN£( 6) 

GRIDS ARE PRESENTLY LIMITED BY DIMENSION SIZE TO 140 POINTS 
IN XI DIRECTION AND 80 POINTS IN ETA DIRECTION. THESE 
DIMENSION SIZES ARE FOUND IN SUBROUTINES INNER, OUTER, 

RELAX, CLUSTR, EPSIL, AND PLAWT AS WELL AS HERE IN THE 
MAIN PROGRAM. 
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c 

DATA PI/3.141592654/ 

C 

C READ IN PARAMETERS. 

RE AD( 5 * 101 ) JWAKE* KMAX* MAXI T 
RE AD( 5*102 ) XGMX,XGMN,YMAX,YMlN,XNOSE, XT AIL 
READ(5*102> DY1 
READ! 5*102 ) XORG* £ T ACD* BET A 
RE AD( 5* 102 ) OMEGA 
R EAO( 5* 102 ) DY2 
CALL RE AO IN 
READ( 5,10GUTITLE 
READ( 5*100 J LINE 
C 

WRITE! 6* 112 ) JWAKE* KMAX»MAXIT* XGMX* XGMN* YMAX* YM1N, XNQSE* XTAIL 

WRITE! 6* 113) 0Y1*X0RG*ET ACD* BETA* OMEGA* DY2 

WRITE!6*137)ITITLE 

WR I TE( 6* 109) LINE 

ETAC«6TACD*PI/130. 

C 

C READ IN AIRFIOL SHAPE. 

N ■ 0 

1 N « N *1 

RE AD! 5* 102 ) XB(N)*YB!N) , 

IF! EOF !5 ) ) 2,1 

2 N BOD » N-l 
WRI TEI 6* 103 ) 

WRITE (6,104) (N,XB (N), YB (N ) , N« 1,N BOD ) 

C 

C DISTRIBUTE POINTS ON INNER AND REAR BOUNDARIES. 

CALL INNER! NBCD, JWAKE* KMAX* XB* YB* XNOSE* XT AIL* XGMX* YMAX* Y MIN, 

1 DYI* JMAX ,X»Y) 

C 

C DISTRIBUTE POINTS ON BOTTOM-FRONT-TOP BOUNDARY. 

CALL OUTER! XGMX, XGMN* YMAX* YMIN, XORG* ETAC, BET A* JM AX, KM AX * X , Y ) 

C 

C FILL-IN BETWEEN ETA-MAX AND ETA-MIN BOUNDARIES. 

34 JMM * JMAX -1 
KMM « KMAX— 1 
RKMM«1 ,0/KMM 
C 

DO 13 J*2, JMM 

DELX-1X! J,KMAX)-X( J,l) )*RKMM 
DELY»(Y(J,KMAX)-Y( J*l) ) *RKMM 
C 

DO 13 K =2* KMM 

X ( J , K) « X(J,1) ♦ ( K-DADELX 
13 Yi J,K ) - Y ! J * 1 ) ♦ i K-l ) *DE LY 
C 

C PLOT INITIAL CONDITIONS. 

CALL PL AWT (JMAX, KMAX, X, Y, XGMX* XGMN* YMAX, YMIN) 

C 

C GENERATE UN-CLUSTERED GRID. 

CALL RELAX (JMAX, KMAX* X,Y, OMEGA, MAXI T) 

C 

C PLOT UN-CLUSTERED GRID. 

CALL PLAWT( JMAX* KMAX, X, Y, XGMX, XGMN* YMAX* YMIN ) 
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c 

C APPLY EXPONENTIAL CLUSTERING. 

CALL CLUSTR! JMAX,X,Y, DY2,KMAX> 

C 

C PLOT EXPONENTIALLY CLUSTERED GRID. 

CALL PLAWT ( JM AX, KMAX, X, Y» XGMX, XGMN, YM AX, YMIN ) 

C 

C WRITE FINAL SOLUTION. 

WRITE (6) ((X( J,K), J«l, JMAX),K-1,KMAX),! !Y!J,K), J-l, JM AX) , K«l, KM AX > 
C 

C PLOT BLOW-UP OF INNERMOST 10 LINES. 

NJW* JWAKE/ 4 
NJM-NBQ0+2+NJW 
JJJ-JWAKE-NJW-1 
XMAX-X! JWAKE, 1) 

XMIN*XMAX ' 

YMAX= Y (JWAKE, 1 ) 

YMI N« YMA X 
C 

DO 35 K*l, 10 
DO 35 J= 1, N JM 
jj. jjj+j 
XX«XI JJ,K) 

YY-Y! JJ,K) 

IF( XX.GT.XMAX) XMAX«XX 
IF(XX.LT.XMIN) XMIN-XX 
IF C YY.GT.YMAX) YM AX«YY 
IFIYY.LT. YMIN) YMIN-YY 
X< J,K)*XX 
35 Y!J,K)»YY 
C 

CALL PLAWT ( NJM,10, X, Y, XN AX, XM IN, YM AX , YMIN) 

CALL EOFTV 
C 

STOP 

C 

IOC FORMAT! 8A1G) 

1C1 FORMAT 1 1615 ) 

102 FORMAT (8F10.0 ) 

103 FORMAT! /27H AIRFOIL POINTS AS READ IN.) 

104 FORMAT ! 3H N«, 13, 5X* 3HXB-, E12 .5, 5X, 3HYB», E12.5 ) 

109 FOR MAT! 34H CARD NO. 9i AIRFIOL DESCRIPTION* , 10X,8A10 ) 

112 FORMAT !24H1 PARAMETERS AS READ IN*// 

1 60H CARD NO. 1* JWAKE, NUMBER OF POINTS IN EACH SIDE OF WAKE «, 

2 I5/14X»41HKMAX, NUMBER OF POINTS IN ETA DIRECTION «,I5/ 

3 14 X, 57HM AX IT » MAXIMUM NUMBER OF ITERATIONS IN MAKING UN-CLUSTERE, 

4 8H0 GRID ■ » I 5/ 

5 57H CARD NO. 2* XGMX, X-DIRECTION COORDINATE OF DOWNSTREAM , 

6 10HB0UNDARY •, F20.1C/14X, 32HXGMN, X-DIRECTION COORDINATE OF , 

7 1 9HUPST RE AM BOUNDARY » , F20 .10 /14 X,25HY MAX, Y-DIRECTION COORDIN, 

8 23HATE OF UPPER BOUNDARY », F 2 0. 10/ 14X , 23HYMI N, Y-DIRECTION COORD, 

9 25HINATE OF LOWER BOUNDARY *, F2 0.10/14 X, 19HXN0SE, X-DIRECTION , 

A 39HC00RDINATE OF LEADING EDGE OF AIRFOIL »# F20. 10/ 14X, 

B 59HXTAIL, X-DIRECTION COORDINATE OF TRAILING EDGE OF AIRFOIL 
C F20.10) 
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113 FORMAT ( 51H CARO NO. 3* DY1, MINIMUM Y-INCP.EMENT ON REARWARD , 

1 33H80UNDARY FOR INITIAL CONDITIONS »#F20.10/1AH CARD NO. A» , 

2 61HXORG, X-DIRECTION LOCATION OF ORIGIN FOR ANGULAR CLUSTERING «» 

3 F20.10/1AX,4J>HETAC, ANGLE (IN DEGREES) ABOUT WHICH ANGULAR , 

A 20HCLUSTER ING IS DONE -t F8. 2/ 1AX, 25HB ETA, ANGULAR CLUSTERING , 

5 11HPARAME1ER «,F23.10/35H CARD NO. 5« OMEGA, PARAMETER FOR , 

6 17HTH0MPS ON SOLVER «,F20.10/31H CARD NO. 6« DY2, MINIMUM ETA-, 

7 AOHDIR ECTION SPACING FOR FINAL CLUSTERING », F20.1Q) 

137 FORMAT (25H CARO NO. 8« PLOT TITLE »# 10X, e AlO ) 

C 

END 
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SUBROUTINE INNER 


Subroutine INNER distributes points on the airfoil-wake (n = n m i n ) 
boundary and the rear boundaries. The distribution of points on the airfoil 
part of the r) = n m i n boundary is taken from the input data cards describing 
the airfoil shape; no redistribution of these points or modifications of their 
number is allowed for in the program. The airfoil shape, however, is shifted 
and scaled to allow arbitrary placement on the a: axis of the leading and 
trailing edges. Points on the wake part of the q = n m i n boundary are exponen- 
tially stretched toward the rearward boundary using the transformation 


S 


1 


= 0 





+ AS. 


mm 



e 

A - 1 



(18) 

(19) 


much like equations (11) and (12) . Distribution of points in the y direction 
on the rear (£ = g mi - n and £ = £ max ) boundaries is also obtained by exponential 
stretching from the wake to the bottom and top boundaries using the transfor- 
mation in equations (11) and (12). It should be noted that the ^-direction 
point distribution on the rear boundaries fa and ei of figure 2 is used only 
as a boundary condition in the TTM method and is discarded along with all 
other q spacing in the application of the final exponential clustering. 


♦ DECK INNER 

SUBROUTINE INNER (NB0D,J WAKE, KMAX,XB,YB, XNOSc, XTAIL, XGMX, YMAX,YM IN, 
1 0Y1»JMAX*X» Y) 

C 

C THIS SUBROUTINE DISTRIBUTES POINTS ON THE ETA-MIN (AIRFOIL- 

C WAKE, INNER) BOUNOARY AND THE REAR (OUTFLOW) BOUNDARY, 

C 

DIMENSION XBC1), YB ( 1 ), X (140, 83 ), Y ( 1A0,80) 

C 

C SHIFT AND SCALE AIRFOIL POIMTS. 

X MI N * XB ( 1 ) 

XMAX « XB ( 1 ) 

DO 3 N*2, NBQD 

I F ( X B(N ) ,LT, X MIN) XMIN • XB(N> 

3 I F ( XB(N) ,GT. XMAX) XMAX ■ X B ( N ) 

C 

SCALE • (XT AI L~X NOSE ) / ( XMAX - XMIN) 

DO A N*l, N BOD 

X B ( N ) • { XB( N) - XMIN)*SCALE*XNOSE 
A YB (N ) » YB(N)*SCAL6 
WRITE! 6, ICS) 

WRITE (6, 106 ) <N,XB(N), YB ( N) , N=l, NBOD > 

C 

JWT « J WAKE ♦ NBOD -1 
N » 0 

00 5 J » JWAKE,JWT 
N ■ N +1 
X(J,1) « XB (N ) 

5 Y ( J , 1 ) - YB ( N ) 
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c 

c 


c 


c 

c 


c 


e 

c 


c 

c 

105 

106 
107 

ioe 

no 

m 

c 


DISTRIBUTE POINTS ON WAKE BOUNDARY, 

DXO ■ ( X ( JWT, 1 ) - X!JWT-1,1> +X ( J WAKE » 1 ) -X ( JW AKE+1, 1 ) ) *. 5 

N « -1 

JMA X « NBOD *2*JWAKE -2 
JWTP • JWT ♦ 1 

XX * 1 X! JWT, 1 ) ♦ X (JWAKE » 1) ) *• 5 
YY • ( Y ( JWT, 1 ) ♦ Y { J WAKE# 1 ) ) * . 5 

EPSW*EPSIL!XGMX,XX,DXO, JM AX-JW T*l, 0 .001# 100, 2, 1 » 

DO 6 J* JWT P, JMAX 
N ■ N +1 

JJ * JWAKE -1 -N 

XX-XX*DX0*11.0+EPSW/SQRT!FLQAT!N+1) > ) **N 

XlJ,l) a XX 

X!JJ,1) • XX 

Y<J,1) * YY 

Y( J J# 1) • YY 

NONE* 1 

WRI TE( 6# 107 ) 

WRITE (6# 108) ( J,NQNE,X< J,l) ,J,NONE,Y( J,l), J«l, JMAX) 

DISTRIBUTE POINTS IN Y DIRECTION ON REARWARD BOUNDARY. 

N«-l 

EP$Y1T«EPSIL(YMAX,YY,OY1»KMAX»0.001»1GO»1»1) 
EPSY1B*EPSIL!-YMIN,-YY,DY1,KMAX,L.001,100,1,2 ) 

DO 8 K=2, KMAX 

N*N*1 

Y( JMAX, Kl-Y! JMAX, K-l ) ♦DY1* < 1 . O+EP S TIT ) *+N 
Y!1,K)«YI1,K-1)-DY1*<1.0*EPSY1B)**N 
X! JMAX,K)-XGMX 
X (1,K )*XGMX 

WRITE (6,110) 

WRITE! 6, 108) (N ONE , K,X ( 1, K ), NONE, K, Y( 1, K) , K-l, KMAX ) 

WRITE <6, 111) 

WRITE! 6, 10 8) ! JMAX,K, X ! JMAX , K ), JMA X, K, Y! JMAX, K ),K» 1, KMAX ) 
RETURN 

FORMAT! /32H AIRFOIL POINTS AFTER RESCALING.) 

FORMAT !3H N« , 13, 5X , 3H XB ■, FI 2 . 5, 5X, 3 HY B«, F12. 5 ) 

FORMAT ! /3QH BOUNDARY ON WAKE AND AIRFOIL*) 

FORMAT! 3H X ( , I 3, 1H, , I 3, 2H ) «, FI 2 . 5, 5X, 2HY( , 13, 1H, , 1 3, 2H) • , F12. 5 ) 
FORMAT ! /33H LOWER PART Of REARWARD BOUNDARY*) 

FORMAT ! /33H UPPER PART OF REARWARD BOUNDARY*) 

END 
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SUBROUTINE OUTER 


Subroutine OUTER distributes points on the bottom-front-top boundary by 
first locating a special "clustering origin" at an arbitrary point on the 
x axis, and then placing the boundary points in an angular fashion about that 
origin. The angular distribution of boundary points is derived from an 
angular clustering transformation given by Davy and Reinhardt (ref. 3). The 
user specifies the angle about which clustering is performed and a parameter 
determining the strength of clustering. Equal angular distribution is avail- 
able as a special case. 


SUBROUTINE OUTER l XM AX, XMIN, Y MAX, YM1N, XORG, ET AC, BET A,4MAX,KMAX,X,Y ) 
C 

C THIS SUBROUTINE PLACES POINTS ON BOTTOM-FRONT-TOP BOUNDARY IN 

C ANGULAR FASHION. 

C 

DIMENSION X(140,80) ,Y(140»80) 

C 

LOGICAL CLUSTR 
C 

DATA PI/3.1 415 92 6 54/ 

C 

SINH(X)*0.5*<EXP(X)-EXP<-X)) 

C 

ETARU«ATAN2(YMAX,XMAX-X0RG) 

ETARL*ATAN2(-YMIN,XMAX-X0RG) 

DETA-(2.0*PI-(ETARU + ETARLn/( JMAX-1) 

CLUSTR * .FALSE. 

I F C BETA. GT. 0.0) CLUSTR - .TRUE. 

IF l • NOT. CLUSTR) GO TO 14 
FACT«PI/<2.0*PI-(ETARU*ETARL) ) 

FACTR-1. 0/FACT 

ET ACT* ( ETAC-E T ARU ) *FAC T 

B*0.5*AL0G( ( 1 •♦( EX P( BE T A) -1, ) *ETACT/ PI ) / < l. + C EXP < -BET A )-l . )* 

1 ETACT/PI)) 

RSB « l./SINH<B) 

14 ETA-ETARU 

ANG1-ATAN2! YMAX, XMIN- XORG) 

ANG2 a ATAN2(YMIN»XMI N— XORG ) *2 .0*PI 
NSIOE « 1 
C 

ETARUD*ETAR U*180 ,/P I 
ETARLD*ETAR L*180. /PI 
ANG1D*ANG1*180./PI 
ANG2D*ANG2*180./PI 

WRITE! 6. 109 ) CLUSTR ,ET ARUD, ETARLD, ANG1D, ANG2D 

IF (CLUSTR. AND. ETAC.LT.ETARU. OR. CLUSTR. AND. ETAC.GT.C 2. 0+PI-E T ARL ) ) 

1 GO TO 22 
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c 

OQ 9 J J*2f JMAX 
J-JMAX+l-JJ 
ETA * ETA ♦ OETA 
IF(.NOT.CLUSTR) GO TO 26 
ETAT«( ETA-ETARU) +FACT 

PHIT-ETACT*(SINH(BETA*E TaT/P I-B)*RS9+1. ) 

PHI -ET ARU+PHIT*F ACT R 
GO TO 27 

26 PHI-ETA 
C 

27 GO TO (It 2>3),NSI0£ 

1 IF(PHI.GT.ANGl) NSIDE-2 

GO TO 3 

2 IF ( PHI • GT »A NG2 ) NSIOE-3 

C 

3 GO TO ( 10f 2 1# 12) f NS ID E 

10 Y(JfKMAX) - Y PAX 

X(J,KMAX» « YMAX /TAN< PHI ) ♦ XORG 
GO TO 21 

11 X(JfKMAX) « XM IN 

Y(J,KMAX) « ( XMIN - XORG) *T AN ( PHI ) 

GO TO 21 

12 X( Jf K M AX ) -X CRG + Y MIN /TAN (PHI) 

Y( J f KMAX) - YM IN 

21 CONTINUE 
C 

ETAD-ETAUEO./PI 

PHI0-PHI+160./PI 

WRITE (6f 113) Jf KMAXf X( Jf KMAX)f Ji KMAXf Y ( Jf KMA X) , E TADf PHIDf NS I Dc 
9 CONTINUE 
C 

RETURN 

C 

22 WR I TE( 6f 11 A ) 

STOP 

C 

109 FORMAT! / 8H CLUST R- » Ilf 5Xf 6HE TARU- f F8 « 2f 5Xf 6HET ARL- , F8 . 2f SXf 

1 5HANGl»f F8.2f5Xf5HANG2-f F8.2//29H OUTER BOUNDARY ON TOPf , 

2 1 8HFR0NTf AND BOTTOM!) 

113 FOR MAT ( 3H X(f I3f lHffI3f2H)-fF12»5f5Xf2HY(fI3flHffI3f2H)-f F12.5f 
l 5Xf AHETA«f Fb,2f5Xf AHPHI-f FS.2f 5Xf 6HNSID6-f III 
11A FORMAT ( /23M ERROR EXIT. BAD ETAC.) 

C 

PND 
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SUBROUTINE RELAX 


Subroutine RELAX creates a grid by solving equations (3) and (4) in the 
uniform £-n domain using successive line overrelaxation (SLOR). Second-order 
accurate central difference operators are used throughout, and tridiagonal 
matrices are solved along lines of constant n using the utility program TRIB. 


♦DECK RELAX 

SUBROUTINE RELAX ( JMAX, KM AX, X# Y, OMEGA# MAXI T) 

THIS SUBROUTINE SOLVES BY SLOR THE DIFFERENTIAL EQUATIONS 
THAT CONSTITUTE THOMP SQN-THAMfcS-MAST IN2S METHOD OF GENERATING 
GRIDS. 

DIMENSION X (140.80) ,Y( 140, 80) 

DIMENSION A{140),B<140> , C <140 ) , D( 140> > F ( 140) , GC140) 

C 

KMM*KM AX—1 
JMM* JMAX— 1 
ICOUNT-O 
C 

2 IC0UNT*IC0UNT+1 
RSUM-0. 

C 

DO I K-2.KMM 
C 

DO 3 J*2, JMM 

XXD»(X(J4l»K)-X(J— 1»K))*0.5 
XED*(X< J,K4l)-X< J,K-1) ) 40.5 
YXD*( Y( J+I,K)-Y( J-1,K) )40.5 
YED«(Y< J,K41)-Y( J.K-lD+0.5 
AD«XED**2+YED**2 
BD«XXD*XED4YXD*YED 
GD"XXD**2+YX0**2 

XX E D* ( X ( J + l.K+l) -X( J + l.K-I >-X < J-1,K41)4X< J-1,K-1) ) +0.25 
YXE0«<Y(j4l,K4l)-Y < J4l,K-l)-Y( J -1, K *1 > 4Y ( J-l, K-i ) ) ♦<). 25 
BD*-2.0*BD 
A( J )« AD 

B( J)*-AO-AD-GD-GD 
C(J )« AD 

F< J>— BD4XXcD-GDMX(J,K*l)4X( J,K-1) ) 

3 G ( J ) *-804YX ED-GD *( Y(J,K4l)4Y<J»K-l) ) 

C 

F (2 )«F(2I-A<2) ♦X(1,K) 

G(2 ) *G (2 )-A ( 2 >*Y ( 1*K) 

F< JMM)*F{ JMM)-C< JMM)4X(JMAX,K) 

G< JMM )*G( JMM )-C< JMM) 4Y( JMAX, K> 

C 

CALL TRIB(A,B,C,0,F,2, JMM) 

CALL TRIB(A,8,C> D, G,2,JMM) 
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c 

DO 4 J*2# JMM 
XC“OMEGA*{F(J)-X(J»K) ) 

YC«OMEGA*<G<J)-Y< J»K)) 

RSUM»RSUM+A8S(XC )+ABS< YC) 

X( J»K)*X< J#K) *XC 
4 Y< J#K)»Y(J#K)+YC 
1 CONTINUE 
C 

WRITE (6» IOC) R SUM# I COUNT 
IF( ICOUNT.LT. MAXIT) GOTO 2 
C 

WRITE (6# 121) 

DO 24 J-1#JMAX>10 
WRITEC6# 115 ) J 

WRITE(6# 11?) ( X ( J#K),K=1»KMAX) 

WR I TE ( 6# 116 ) J 

24 WRITE (6# 117) ( Y ( J# K ) # K«l, KM AX ) 

C 

RETURN 

C 

100 FORMAT C 20H SUM OF RESIDUALS ■ »F20.10# 

1 7H AFTER #I5#12H ITERATIONS.) 

115 FORMAT < /31H X>S FOR CONSTANT XI LINS AT J ■# 15) 

116 F0RMAT(/31H Y>S FOR CONSTANT XI LINE AT J*#I5> 

117 F0RMAT(10E13.6) 

121 FORM AT( /48H AFTER THOMP SON-SOLVER# BEFORE FINAL CLUSTERING!) 

END 
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SUBROUTINE CLUSTR 


Subroutine CLUSTR applies along each line of constant £> the exponential 
clustering transformation discussed at length in the main body of this paper. 


♦DECK CLUSTR 

SUBROUTINE CLUSTR <J MAX, X#Y,DY 2» KMAX) 

THIS SUBROUTINE APPLIES TNE EXPONENTIAL CLUSTERING 
TO THE LINES OF CONSTANT XI. 

DIMENSION X(14C,80),Y(140,80) 

DIMENSION T<80), S ( 80 ) ,XXX < 80) , YY Y < 83) 

C 

KM AX - KMAX-1 
KMM « KMAX-1 
C 

00 20 J-l.JMAX 
C 

Til) * 0. 

DO 16 K*2» KMAX 

16 T(K) ■ TIK-l ) ♦ SORT ( < XU,K) -X ( J, K-l ) ) **2 +( Y( J,K)-Y( J»K-1) ) 
1 * + 2 ) 

C 

EPS Y2-EPSIL(T(KMAX), 0.0, 0Y2, KMAX, 0.001,100, l.J) 

C 

S(1)*0.0 

N« — 1 

DO 15 K c 2» KMAX 
N«N+1 

15 S(K)*S(K— 1)+DY2*(1.0*EPSY2)**N 

C 

DO 17 K«1 » KMAX 
XXX(K) - X<J,K) 

17 YYY(K) * Y(J,K) 

C 

PTL ' 0, 

00 18 K»2,KMM 

18 CALL TAINT(T»XXX,S(K),XU»K),KMAX,2,NER,PTL> 

C 

PTL « 0. 

DO 19 K*2,KMM 

19 CALL TAINT(T,YYY»S(K),Y(J,K),KMAX,2»NER»PTL) 

C 

20 CONTINUE 
C 

WRITE<6,120) 

DO 23 J«1,JMAX 
WRITE ( fc, 115) J 

WR I TE 1 6# 117 ) (X (J,K),K«1,KNAX) 

WRITE(6,1U> J 

23 WR ITE (6/ 117 ) ( Y ( J , K ) , K«1 , KM A X ) 

RETURN 
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c 

120 FORMAT {/2AH AFTER FINAL CLUSTERING*) 

115 FORMAT! / 31H X>S FOR CONSTANT XI LINE AT J«»I5> 

Ufc FORMAT ( /31H Y>S FOR CONSTANT XI LINE AT J«>15) 

117 FORMAT ( 10E 13 • 6 ) 

C 

END 
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FUNCTION EPSIL 


Function EPSIL applies a Newton-Raphson iterative technique to find a 
value of e for use in the exponential stretching transformations given by 
equations (11), (12), (18), and (19). 

FUNCTION £PSIL(FNX,fMIN,DFM,NPT,FPCC,ICC,KEY,NCALL> 

C 

C THIS SUBROUTINE APPLIES A NEVTON-RAPHSON ROOT-FINDING 

C TECHNIQUE TO FIND A VALUE OF EPSILON FOR A PARTICULAR USE 

C OF THE EXPONENTIAL STRECHING TRANSFORMATION. 

C 

DIMENSION R ( 140 ) 

C 

FMXL-FMX 

FMINl-FNIN 

dfml-dfm 

FPCCL-FPCC 

ICCL-ICC 

C 

GO TO <1,2»,KEY 
C 

1 FNPTM2-NPT-2 

IF(NCALL.E0.1 ) EPS»(FHXL/DFHL>**(1.0/FNPTN2)-1.C 
C 

DO 3 NIT-1, ICCL 
EPl-EPSn.U 
EPiTN«EPl**FNPTN2 
REPS-1. 0/EPS 
DFM06-DFML*REPS 

F-FMXL-FMINL-DFNOE*(EP1TN*cP1-1.0> 

IF(A0S(F).LT.FPCCL) GO TO 4 
CFM0c2-DFN0E*REPS 

FPN-0FNOE2*(1.0*EPlTN*(EPS*FNPTM2-1.0) ) 
eps«eps*f/fpn 

3 CONTINUE 

GO TO 5 
C 

2 NPTM-NPT-1 
IF(NCALL.EQ.l J 

1 EPS* < (FMXL/DFMU**C1.0/(NPT-2I )-1.0)*SQRT( FLOATINPTH) ) 

DO 6 L ■ 1, NPTH 

6 R(L)»l.C’/SORT(FLOAT(L) ) 

C 

DU 7 NIT-1, ICCL 

SUM1-C.0 

SUM2-C.0 

DU 8 L-l, NPTM 

FLH2-L-2 

FACTi»l.L*EPS*R(L> 

FACT2-FACT1**FLM2 
SUMl»SUNH-FACT2*FACTi 
8 SUH2-SUN2+(L-D*FACT2*R(L) 

F»FMXL-FMINL-DFNL*SUN1 
IF(A8S(F).LT.FPCCL) GO TO 4 
FPN«0FNL*SUK2 

eps-eps+f/fpn 

7 CONTINUE 

c 

5 EPSIL-EPS 

WRITE (6,1 3D) 

RETURN 

C 
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C 

100 

101 

c 


EPSIl-EPS 

WRIT£(6sl)l) EP$IL,F,NIT 
RETURN 

FORMAT ( /42H EXCEEDED MAX. NO. OF ITERATIONS IN EPSIL.) 
FORMAT ( /7H EPSIL«» F12.5.5X.7H AND F«,F12.5» 5X* 7H AFTER »I3* 
♦ 12H ITERATIONS.) 

END 
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SUBROUTINES TAINT AND TRIB 


Subroutine TAINT is a utility subroutine that inputs a tabulated function 
of one independent variable, then outputs the value of the function at a given 
value of the independent variable by passing a polynomial of arbitrary order 
(<9) through the nearest input points. It is used for second-order 
interpolation. 

Subroutine TRIB solves a tridiagonal system of linear equations using the 
so-called Thomas algorithm (i.e., Gaussian elimination or L-U decomposition). 


SUBROUTINE TAINT ( XTAB, FT AB, X, FX, N,K ,NER » MON) 

SYSTEM LIBRARY SUBROUTINE TAINT FOR POLYNOMIAL INTERPOLATION 
OF A TABULATED FUNCTION. 

DIMENSION XTAB(1),FTAB(1),T< 10>,C(10> 

CPS 0400 TAINT SUBROUTINE- IN FORTRAN II. 

IF (N - K) 1 , 1,2 

1 NER* 2 
RETURN 

2 IF (K-9) 3, 3,1 

3 IF ( MON) A, A, 5 

5 IF ( MON-2) 6,7, A 
A J*0 

NMl-N-1 
DO 6 1*1. NM1 

IF <XTA8U)— XTABU + l) ) 9.11.10 

11 NER-3 
RETURN 

9 J-J-l 
GO TO 6 
10 J-J+l 
8 CONTINUE 
M0N*1 

IF <J) 12.6.6 

12 M0N-2 

7 DC 13 I-l.N 

IF (X-XTAB(I)) 1A.1A.13 
1A J-I 

GO TO 18 

13 CONTINUE 
GO TO 15 

6 DO 16 I *1. N 

IF (X-XTAB(D) 16.17,17 

17 J«I 

GO TO 18 
16 CONTINUE 
15 J*N 

18 J«J— {K+l)/2 

IF (J) 19,19,20 

19 J-l 

20 M* J*K 

IF (M-N) 21,21,22 
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22 J-J-l 

GO TO 2G 
21 KP1-K+1 
JSAVH-J 

26 00 23 L*l, KP1 
C<L)-X-XTAB(J> 

T( L ) ® FT A8 ( J ) 

23 J-J+l 

DO 24 J«1,K 
I»J*1 

25 T(I)=(C(J)*T(I>-CCI)*T( J) >/(C( JI-C(I) ) 
I*I*1 

IF (I-KPl) 25,25 ,24 

24 CONTINUE 


FX«T(KP1) 

NER-1 

RETURN 

END 


SUBROUTINE TRI8( A,8,C, X, F, NL, NU) 

DIMENSION A (2 ),B(2) ,C(2),X<2),FC2) 

THIS SUBROUTINE SOLVES A TRI-DIAGONAL SYSTEM OF LINEAR 
EQUATIONS. 

XCNL)«CCNL)/B(NL» 

F(NL»-F(NL)/B(NL) 

NLP1 « NL +1 
DO I J«NLP1»NU 
Z«1«/(R( J)-A( J)+X( J-m 
X«JI»C(J)*Z 

1 F ( J > * « F ( J ) -AC J)*F(J-1))*Z 
NUPNL-NU+NL 

DO 2 Jl a NL PI, NU 
J-NUPNL-J1 

2 F( J)«F( J»-X(J )*F (J+l) 

RETURN 

END 
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SC4020 PLOT SUBROUTINES 


Subroutines PLAWT, AXIS, LABEL, PLOT, and TITLE are used to illustrate 
the grids on the SC4020 plotter. They contain calls to the standard SC4020 
software subroutines; thus, the SC4020 software subroutines must be made avail- 
able to the loader. If the user's computer installation does not include an 
SC4020, the five subroutines listed above should be removed from the deck; 
calls to subroutines READ IN, PLAWT, and EOFTV should be removed from the main 
program; and data card 7 should be deleted. 


SUBROUTINE PLArfT(N,H,X,Y,XMAX#XMIN, YMAX.YMIN) 

C 

C THIS SUBROUTINE IS ONE OF THE SUBROUTINES THAT PLOTS GRIDS 

C ON THE SC4O20 PLOTTER. 

C 

DIMENSION X(14O»8O)»Y(14O»0O)»XX(8O),YY(8O) 
COMMON/TCOM/ITITLE(fl) 

C 

C READJUST PLOT LIMITS SO AS TO AVIOD A STRECHED PLOT. 

XOIF«XMAX-XMIN 
YDIF-YMAX-YMIN 
IF C XDIF.LT.YDIF) GO TO 4 
XDIFH«XDIF*0.5 
YMID=(YMAX4YMIN)*0.5 
YMX-YMID4X0IFH 
YMN *YMI D-XOIFH 
XMX »X MAX 
XMN-XMIN 
GO TO 5 

4 YDI FH a Y D IF*0. 5 

XMID«< XMAX+XMIN) *3.5 
XMX*XMlD+YDIFri 
XMN-XMID-YDIFH 
YMX-YMAX 
YMN-YMIN 
b CONTINUE 

C 

C PLOT THE LINES. 

CALL AXIS<XMN,XMX,YMN, YMX,0,0) 

CALL TI TLE ( ITI TLE, 80 ) 

001 J-1,M 

1 CALL PL0T( X(l, J>, Y(1,J ),N,0) 

D02 1*1* N 

D03 J*l» M 
XX< J)"XCI, J) 

3 YY ( J ) *Y ( I, J ) 

2 CALL PLOTC XX, YY»M,0> 

RETURN 

END 
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SUBROUTINE AXIS( XMIN# XM AX# YM IN# YM AX # NXI NTS » NY I NTS ) 

THIS SUBROUTINE IS ONE OF THOSE THAT PLOTS GRIDS 
ON THE SC4020 PLOTTER. 

THIS SUBROUTINE ORAWS THE X ANO r AXES# MARKS OFF THE 
INTERVALS# SETS THE SCALE VALUES# AND MOVES THE FILM TO A NEW 
FRAME. 

XM I N IS THE MINIMUM VALUE OF THE X (HORIZONTAL) VARIABLE. 

XMAX IS THE MAXIMUM VALUE OF THE X VARIABLE. 

YMIN IS THE MINIMUM VALUE OF THE Y (VERTICAL) VARIABLE. 

YM AX IS THE MAXIMUM VALUE OF THE Y VARIABLE. 

NXINTS IS THE NUMBER OF INTERVALS INTO WHICH THE X AXIS IS TO 
BE DIVIDEO. 

NYI NTS IS THE NUMBER OF INTERVALS INTO WHICH THE Y AXIS IS TO 
BE DIVIDED. 

DO INITIALIZATION 
CALL C AMRAV( 35) 

CALL SMALL V 
CALL FRAMEV(O) 

NX-NXINTS 
NY* NYINTS 

C IF NX OR NY • LE. 0 MAKE THEM EQUAL TD 1 

IF (NX .LE. 0) NX-1 
IF (NY .LE. 0) NY-1 

C SET UP VALUES TO BE USEO FOR MARGINS 
ML-123 
MR-923 
MB- 123 
MT-923 

C DRAW X AXIS 

CALL LINEV(ML#MB,MR,MB) 

C DRAW Y AXIS 

CALL LINEV(ML#MB#ML#MT) 

C DETERMINE INCREMENTS FOR TIC MARKS 
DX«( XMAX-XMINJ/NX 
DY- (YMAX-YMINI/NY 
C SCALE X ANO Y VALUES 

CALL XSCALV(XMIN,XMAX#ML»100) 

CALL YSCALV(YMIN#YMAX#M8#100> 

C DRAW TIC MARKS ON THE X AXIS 

CALL LINRV(1#M8-20#MB-5#M8 + 5#XMIN#XMAX# DX# 0# -1# -3# 8 ) 

C DRAW TIC MARKS ON Y AXIS 

CALL LINRV(2»ML- 9D» ML-5# ML+ 5# YMIN# YMAX# DY# 0# -1 #-3# 10 ) 

RETURN 

END 
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SUBROUTINE LABEL < IXLABL,NX, I YLABL, NY) 

THIS SUBROUTINE IS ONE Of THOSE THAT PLOTS GRIDS 
ON THE SC4C20 PLOTTER. 

THE MAXIMUM NUMBER OF CHARACTERS ALLOWED FOR THE X LABEL IS 54. 
THE MAXIMUM NUMBER OF CHARACTERS ALLOWED FOR THE Y LABEL IS 36. 
IN X* IAB S (NX ) 

INY* I AB S( NY ) 

IF (INX .GT. 54) INX*54 
IF (INY .GT. 38) INY-38 
IX* 500-( I NX /2 ) *1 8 
IY*51Q+(lNY/2)*26 
LABEL THE X AXIS 

CALL RITE2V( IX. 60. 1000. 90# 1. INX, 1, IXL ABL. NL AST ) 

LABEL THE Y AXIS 
DO 100 1*1, INY 

ITEMP«IY-26*<I-1) 

CALL R ITE2V( 10, I TEMP, 50, 90,1, 1, I, 1YL ABL, NLAS T) 

ICC CONTINUE 
RETURN 
END 


SUBROUTINE PLOT( X , Y, NBR » NSYM ) 

THIS SUBROUTINE IS ONE OF THOSE THAT PLOTS GRIDS 
ON THE SC4020 PLOTTER. 

DIMENSION X(1)»Y(1),MARKPT(5) 

DATA MARKPT / 42, 16, 55, 38,44/ 

C SYMBOLS ARE, IN ORDER* , + X 0 * 

IF (NSYM .GT. 0 .AND. NSYM .LE. 5) GO TO 100 
J* I ABS ( NBR) -1 
DO 110 1*1, J 

CALL L INE V( IXV(X(I>), IYV(Y(I)>, I XV( X ( I+l ) ) , IYV ( Y ( I+l) ) ) 
11C CONTINUE 
RETURN 

100 CALL APLOTVtNBR, X, Y, 1,1,1, MARKPT (NSYM ) ) 

RETURN 

END 
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SUBROUTINE TITLE (ITITLE.NCHARS) 


THIS SUBROUTINE IS ONE OF THOSE THAT PLOTS GRIDS 
ON THE SC402D PLOTTER. 

THE MAXIMUM NUMBER OF CHARACTERS ALLOWED IN THE TITLE IS 108. 

ICHARS«IABS(NCHARS) 

IF CICHARS .GT. 108) ICHARS-106 
IF (ICHARS .GT. 54) IX = 14 
IF (ICHARS .LE. 54) IX • 510-< I CHAR S tl ) *1 8 
I Y* 990 

CALL RITE2V(IX.IY»101C»9C.l, ICHARS>1» ITITLE. NLAST) 

RETURN 

END 



INPUT DATA CARDS 


Input data cards are as follows. All integers are read using 1615 format. 
All floating-point variables are read using 8F10.0 format. All alphanumeric 
data are read using 8A10 format. 


Card 1. 

Card 

1 contains three integers. 



Columns 

Name 

Description 



1 to 5 

JWAKE 

Number of grid points in the x direction in 
including the trailing edge point 

the 

wake 

6 to 10 

KMAX 

Number cf grid points in the r> direction 



11 to 15 

MAXET 

Number of iterations to be used in creating 
tered grid 

the 

unclus- 


Card 2 . Card 2 contains floating point numbers that define the physical 
size of the grid and the location of the airfoil. 

Columns Name Description 

I to 10 XGMX x-direction coordinate of rearward boundary 

II to 20 XGMN ^-direction coordinate of front boundary 

21 to 30 YMAX ^-direction coordinate of top boundary 

31 to 40 YMIN ^-direction coordinate of bottom boundary 

41 to 50 XNOSE tf-direction coordinate of leading edge of airfoil 

51 to 60 XTAIL x-direction coordinate of trailing edge of airfoil 

Card 3 . Card 3 contains in columns 1 to 10 one floating-point number 
DY 1, which is the physical spacing in the y direction on the rear boundary 
immediately above and below the wake. 

Card 4 . Card 4 contains floating-point numbers determining the location 
of points on the bottom-front-top boundary. 

Columns Name Description 

I to 10 XORG Location on the x axis of the special "clustering origin" 

II to 20 ETACD Angle (in degrees), measured counterclockwise from the 

X axis, about which clustering should be performed 

21 to 30 BETA Parameter determining the strength of clustering. Zero 

(for equal-angular distribution) or values in the range 
1 to 5 are recommended. 

Card 5 . Card 5 contains one floating-point number in columns 1 to 10, 
OMEGA, the relaxation parameter for the SLOR procedure in subroutine RELAX. 
Values must be between 0 and 2, with a typical value being 1.55. 
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Card 6 . Card 6 contains one floating-point number in columns 1 to 10, 
the parameter DY 2, which is the physical spacing in the p direction between 
the airfoil-wake boundary and the next line of constant p. 

Card 7 . Card 7, read by SC4020 subroutine READ IN, contains up to 80 
alphanumeric characters to be used for a comment on the title frame of the 
SC4020 plots. 

Card 8 . Card 8 contains up to 71 alphanumeric characters to be placed as 
a title on each of the four SC4020 plots that result from each run. 

Card 9 . Card 9 contains up to 80 alphanumeric characters describing the 
airfoil to be used. This description will appear on the printed output. 

Cards 10 , 11 , . . . . All cards following card 9 contain in columns 1 to 10 
and 11 to 20 two floating-point numbers, x and y, respectively, which are 
coordinates of points on the airfoil surface, defining the airfoil shape. The 
points should be ordered to begin and end at the trailing edge and proceed 
counterclockwise . 

Two notes of caution with respect to the input data cards are in order 
for those installing this program on computers of manufacture other than 
Control Data. First, the cards defining the airfoil shape are read until an 
end-of-file is encountered. End-of-file conditions are handled differently by 
different computers. Second, alphanumeric data is stored 10 characters per 
word on Control Data machines. Care should be taken on machines having dif- 
ferent word lengths. 


SAMPLE DATA CASE 


The data cards used to produce the grid shown in figure 5 follow: 


12 28 

2C0 

8.0 

6.0 

0.04 


1.55 


0.01 


REESE SORENSON NASA 

NACA 0012 


NACA 0012 


1.000000 

-.001251 

.981540 

-.003816 

.959694 

-.006776 

.929861 

-.010690 

.892711 

-.015378 

.848960 

-.020639 

.799449 

-.026266 

.745313 

-.032084 

.687413 

-.037853 

.626854 

-.043365 


— 6*0 


STOP 5124 


0.0 1.0 


PH. 9b 5-5125 TAPE CO 
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•564684 
.501943 
.439612 
.378638 
.319904 
.264227 
.212357 
.164970 
.122669 
.085973 
.055327 
.031109 
.013637 
.003244 
-.000000 
.003283 
.013651 
.031111 
.055331 
.0 85 975 
.122670 
.164972 
.212357 
.264227 
.319904 
.378638 
.434613 
.501943 
.564685 
.626655 
.687413 
.745314 
.799450 
.848961 
.892711 
.929881 
.959694 
.981540 
1.00CC00 


-.048424 
-.052 817 
-.056319 
-.058746 
-.059927 
-.059712 
-.056027 
-.054863 
-.050245 
-.044279 
-.037107 
-.028668 
-.019743 
-.009929 
.000002 
.009932 
.019726 
.028866 
.037100 
,044274 
.050244 
.054 656 
.058027 
.054712 
.059926 
.058746 
.056319 
.052815 
.048424 
.043366 
.037653 
.032064 
.026269 
.020640 
.015381 
.010691 
.006777 
.003816 
.001252 
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Figure 2.— Schematic 
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(a) Overview. 

Figure 3.— Computed grid using equations (3) and (4) with specified 

boundaries . 
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(b) Grid detail near boundary. 
Figure 3.— Concluded. 
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Figure 4.— Detail of computational grid using equations (7) and (8). 
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(b) Grid detail near boundary. 

Figure 5.— Computational grid after application of clustering along 

lines of constant 
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